# Imports
from coopr.pyomo import *
from coopr.opt import SolverFactory
from ReferenceModel import model
import numpy

# Solve EV for given number of sample realizations with fixed X at X_EV
numSamples = 500
numX=5
optVal=numpy.array([0 for i in range(numSamples)])

# Choose the solver
opt = SolverFactory('gurobi')

# See the result from part b
EV_X = [9.2, 23.45, 5.1, 8.3, 4.95]

# Iterate through all samples
for i in range(numSamples):
    datafile = './scenariodata/Scenario' + str(i+1) + '.dat'
    instance = model.create(datafile)

    # Fix values of x at x_EV 
    for j in range(numX):
        instance.X[j+1] = EV_X[j]
        instance.X[j+1].fixed = True
    instance.preprocess()

    # Solve the instance 
    results = opt.solve(instance)
    print "Solve" + str(i) + "th instance"
    #print instance.X.extract_values()
    instance.load(results)
    optVal[i] = value(instance.TotalProfit)
    
# Point estimate
EEV = optVal[:].mean()

# Interval estimate
z_alpha = 1.96
EEV_var = optVal[:].var()*numSamples/(numSamples-1)
EEV_halfwidth = z_alpha * sqrt(EEV_var/numSamples)
EEV_CI_lo = EEV - EEV_halfwidth
EEV_CI_hi = EEV + EEV_halfwidth
    
print EEV
print EEV_CI_lo, EEV_CI_hi

